%% we define function to find fixed points
function [] = findFixedPoints(index, nLS, nNLS, nMB, nFB,rhob,rhoo)
    global alpha qn qb eta kappa psim psiw r_initial  v0 s_output markup_output add_mp_output std_mp_output %outside option: v0
  
    nplayers = nLS+nNLS+nMB+nFB;
    tol = 10^(-5);
    MaxIterations =1000;
    iters = zeros(nplayers,1);
    diffs = 1000*ones(nplayers,1);

    % create initialized values only for this player
    if index > nplayers
        return;
    end
    
    
    while diffs(index)>tol && iters(index) < MaxIterations
        findFixedPoints(index+1, nLS, nNLS, nMB, nFB,rhob,rhoo)

        % do one iteration of calculation
        iters(index) = iters(index) + 1;
        
        % shadow bank financed by full function banks
        Als = exp(-alpha*r_initial(1:nLS) +qn);%exp((-alpha*r_initial(1:nLS) +qn)/lambdaF) -- nested logit
        Zls  = sum(Als);

        % shadow bank financed by warehouse banks
        Anls = exp(-alpha*r_initial(nLS+1:nLS+nNLS) +qn);%exp((-alpha*r_initial(nLS+1:nLS+nNLS) +qn)/lambdaF) -- nested logit
        Znls  = sum(Anls);
        
        % mortgage bank
        Amb = exp(-alpha*r_initial(nLS+nNLS+1:nLS+nNLS+nMB) +qb);%exp((-alpha*r_initial(nLS+nNLS+1:nLS+nNLS+nMB) +qb)/lambdaF) -- nested logit
        Zmb  = sum(Amb);     
        
        % full function bank
        Afb = exp(-alpha*r_initial(nLS+nNLS+nMB+1:end) +qb);%exp((-alpha*r_initial(nLS+nNLS+nMB+1:end) +qb)/lambdaF) -- nested logit
        Zfb  = sum(Afb);     
        
        Dem = v0+Zls+Znls+Zmb+Zfb;

        sLS = Als ./ Dem;
        sNLS = Anls ./ Dem;
        sMB = Amb ./ Dem;
        sFB = Afb ./ Dem;

        r_LS = eta*rhob+kappa+1./(alpha*(1-sLS));
        r_NLS = eta*rhoo+kappa+1./(alpha*(1-sNLS));
        r_MB = psim+1./(alpha*(1-sMB));
        r_FB = psim+1./(alpha*(1-sFB))+eta.*(rhob-psiw).*sum(sLS)./(1-sFB); 
        
        r = [r_LS;r_NLS;r_MB;r_FB];
        s = [sLS;sNLS;sMB;sFB];
      
        std_markup = [1./(alpha*(1-sLS));1./(alpha*(1-sNLS));1./(alpha*(1-sMB));1./(alpha*(1-sFB))];
        add_markup = [zeros(nLS,1);zeros(nNLS,1);zeros(nMB,1); eta.*(rhob-psiw).*sum(sLS)./(1-sFB)];

        markup = [1./(alpha*(1-sLS));1./(alpha*(1-sNLS));1./(alpha*(1-sMB));1./(alpha*(1-sFB))+eta.*(rhob-psiw).*sum(sLS)./(1-sFB)];
      
        diffs(index) = abs(r(index)-r_initial(index));
        r_initial(index) = r(index);
        s_output(index) = s(index);
        markup_output(index) = markup(index);
        add_mp_output(index) = add_markup(index);
        std_mp_output(index) = std_markup(index);
        
       % display(diffs);


            
    end
end